disp('flf');
an = @(x, n) 1/2/pi*n./((n/2).^2 + x.^2);

x = 1.2;
kc = 2;

format long;
disp(fn_f2(x, kc));
disp(fn_f4(x, kc));

n = rand();
kc = rand();
k = rand();

disp(fn_int_tail(k, kc, n));
disp(fn_f2(-2*k/n, 2*kc/n)*8/2/pi/n^2);
disp(fn_f2(-2*k/n, 2*kc/n)*8/2/pi/n^2 - fn_int_tail(k,kc,n))

v2 = rand(1,10);
v4 = rand(1,10);
v = rand(1,10);

[alpha, beta] = fn_fit(v2, v4, v);
tv = v - alpha*v2 - beta*v4;
fprintf('alpha = %f  beta = %f \n', alpha, beta);
disp(tv*transpose(v2));
disp(tv*transpose(v4));

vx = linspace(0, 10, 400);
v_an2 = an(vx, 2);
v_an3 = an(vx, 3);

v_an5 = an(vx, 5);
vx_out = vx;

[vf_out0, vf_out2, vf_out4] = fn_anf(3, vx, v_an2, vx_out);

close(figure(1));
figure(1);
plot(fn_ex1(vx), fn_ex2(v_an5), '-r');
hold on;
plot(fn_ex1(vx), fn_ex2(vf_out0), 'o');
plot(fn_ex1(vx), fn_ex2(vf_out2), '*');
plot(fn_ex1(vx), fn_ex2(vf_out4), 's');

vx = linspace(0,10,1000);
vx_half = vx(1:500);
vx_half2 = vx(501:1000);
vf = an(vx, 2);
vf_half = an(vx_half, 2);
vf_half2 = an(vx_half2, 2);


tf = vf_half;
v_tail = vx_half(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
tf_tail = vf_half(end-9:end);
[alpha, beta] = fn_fit(v2, v4, tf_tail);

close(figure(2));
figure(2);
plot(vx_half2, vf_half2, '-r');
hold on;
plot(vx_half2, alpha./vx_half2.^2, 'o');
plot(vx_half2, alpha./vx_half2.^2 + beta./vx_half2.^4, 's');

alpha0 = vx_half(end)^2*vf_half(end);
plot(vx_half2, alpha0./vx_half2.^2, '*');
fprintf('alpha0 = %f \n', alpha0);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function res = fn_yn_inf(n, h0)

if h0>0
    
    res = n*h0 ...
        + log((1+exp(-(2*n+2)*h0)-exp(-n*h0)-exp(-(n+2)*h0)) ...
        ./(1+exp(-2*h0)-2*exp(-h0)));

else
   
    res = log(n.^2 + 2*n);
    
end

end


function vu_out = fn_Tnm(n, m, a, vx, vu, vx_out)

vu_out = 0;

vf = log((1 + exp(-vu)*a)./(1 + a));

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
        
    end
    
end

end




function res = fn_inveta(m, h0)

if h0>0
    res = sinh(h0/2)^2/(sinh(h0/2*(m+1))^2 - sinh(h0/2)^2);
else
    res = 1/(m^2 + 2*m);
end
    
    
end



function [vf_out0, vf_out2, vf_out4] = fn_anf(n, vx, vf, vx_out)

an = @(x) 1/2/pi*n./((n/2).^2 + x.^2);

vk = vx_out;
vp = vx;

mk = transpose(vk)*ones(1,length(vp));
mp = ones(length(vk),1)*vp;
ma = an(mk - mp) + an(mk + mp);

vpd = vp(2:end) - vp(1:end-1);
vp_nt = zeros(1,length(vp));
vp_nt(1:end-1) = 1/2*vpd;
vp_nt(2:end) = vp_nt(2:end) + 1/2*vpd;

vf_out = ma*diag(vp_nt)*transpose(vf);
vf_out = transpose(vf_out);
vf_out0 = vf_out;

kc = vp(end);
v_tail = vp(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

fprintf('alpha = %f, beta = %f \n', alpha, beta);

vf_out_addon2 = alpha ...
    *8/2/pi/n^2*(fn_f2(2*vk/n, 2*kc/n) + fn_f2(-2*vk/n, 2*kc/n));

vf_out_addon4 = beta ...
    *32/2/pi/n^4*(fn_f4(2*vk/n, 2*kc/n) + fn_f4(-2*vk/n, 2*kc/n));

vf_out2 = vf_out0 + vf_out_addon2;

vf_out4 = vf_out0 + vf_out_addon2 + vf_out_addon4;

end


function vx_ex1 = fn_ex1(vx)

vx_ex1 = [-vx(end:-1:2), vx];

end

function vx_ex2 = fn_ex2(vx)

vx_ex2 = [vx(end:-1:2), vx];

end


function res = fn_int_tail(k, kc, n)

res = 8/pi*n*log(4)*k./(4*k.^2 + n^2).^2 ...
    + (4*k.^2 - n^2)*2./(4*k.^2 + n^2).^2 ...
    + 2/pi*n./(4*k.^2 + n^2)./kc ...
    - 4/pi*(4*k.^2 - n^2)./(4*k.^2 + n^2).^2.*atan(2*(k + kc)./n) ...
    + 8/pi*n*k./(4*k.^2 + n^2).^2.*log(kc.^2./(n^2 + 4*(k + kc).^2));

end


function res = fn_f2(x, kc)

res = -pi/2*(1 - x.^2)./(1 + x.^2).^2 ...
    + 1/kc./(1 + x.^2) + (1 - x.^2)./(1 + x.^2).^2.*atan(kc - x) ...
    + x./(1 + x.^2).^2.*log((1 + (kc - x).^2)./kc^2);

end


function res = fn_f4(x, kc)

res = 3*pi*(1 - 6*x.^2 + x.^4) ...
    + 6*(1 - 6*x.^2 + x.^4).*atan(x - kc) ...
    + 12*x.*(x.^2 - 1).*log((1 + (x - kc).^2)./kc^2) ...
    + 2*(1 + x.^2).^3/kc^3 ...
    + 6*x.*(1 + x.^2).^2/kc^2 ...
    + 6*(1 + x.^2).*(3*x.^2 - 1)/kc;

res = res./6./(1 + x.^2).^4;

end


function [alpha, beta] = fn_fit(v2, v4, v)

m = [v2*transpose(v2), v2*transpose(v4); ...
    v4*transpose(v2), v4*transpose(v4)];

RHS = [v2*transpose(v); v4*transpose(v)];

tv = m\RHS;
alpha = tv(1);
beta = tv(2);

end
